home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / examples / estimate.r < prev    next >
Text File  |  1994-04-25  |  4KB  |  153 lines
  1. // -------------------------------------------------------------------
  2. //
  3. // Beginning of codes to test memory allocation/deallocation.
  4. //
  5. // The following lines are in a file called `estimate.r'
  6. initstate = function(v) 
  7. {
  8.   local (prob, tmp, i);
  9.  
  10.   tmp = 0;
  11.   prob = rand();
  12.   for(i in 1:v.nr) {
  13.     tmp = tmp + v[i;1];
  14.     if (tmp >= prob) {
  15.       break;
  16.     }
  17.   }
  18.   if (v[i;1] == 0) {
  19.     error("initstate: Attempting to initialize to invalid state.");
  20.   }
  21.   return i;
  22. }
  23.  
  24. nextstate = function(thisstate, trans) 
  25. {
  26.   local (prob, tmp, i);
  27.  
  28.   tmp = 0;
  29.   prob = rand();
  30.   for(i in 1:trans.nc) {
  31.     tmp = tmp + trans[thisstate;i];
  32.     if (tmp >= prob) {
  33.       break;
  34.     }
  35.   }
  36.   if (trans[thisstate;i] == 0) {
  37.     error("nextstate: Attempting transition to invalid state.");
  38.   }
  39.   return i;
  40. }
  41.  
  42. randwalk = function(A, P, f, init, h, p, chlen) 
  43. {
  44. //
  45. // Randwalk uses the transition matrix P, the initial state init,
  46. // and the number of steps chlen to take a random walk.  The other
  47. // input variables are required in order to keep the statistics we
  48. // need to accumulate. 
  49. //
  50.   local (old, new, i, W, total);
  51.  
  52.   old = init;
  53.   W = 1.0;
  54.   total = f[old;1];
  55.   for (i in 1:chlen) {
  56.     new=nextstate(old, P);
  57.     W = W * A[old;new] ./ P[old;new];
  58.     total = total + W * f[new;1];
  59.     old = new;
  60.   }
  61.   return h[init;1] * total ./ p[init;1];
  62. }
  63.  
  64. estimate = function(B, f, h, chlen, niters) 
  65. {
  66. // Given a matrix equation Bx = f with known coefficient matrix B and
  67. // known right-hand-side f, the function `estimate' uses a Markov chain
  68. // technique to approximate the inner product of the known vector h
  69. // with the unknown vector x (or rather, the inner product of h with the
  70. // chlen+1st iterate of the fixed point iteration x = Ax+f, where
  71. // B = I - A).  The input `niters' controls how many random walks are taken.
  72. // For details of the method, see Chapter 5 of R. Rubinstein's ``Simulation
  73. // and the Monte Carlo Method.''
  74.   local (A, n, i, j, tmp, P, p, total, eta);
  75.   n = B.nr;
  76.   A = eye(n,n) - B;
  77. //
  78. // Now generate P, the ergodic Markov chain, by scaling the rows of A so
  79. // that the row-sums are all = 1.  Note: we are making an implicit assumption
  80. // that the infinity norm of the matrix A is less than one.  (see Rubinstein).
  81. //
  82.   P=zeros(n,n);
  83.   for(i in 1:n) {
  84.     tmp = mysum(abs(A[i;]));
  85.     P[i;] = abs(A[i;])./tmp;
  86.   }
  87. //
  88. // Generate the initial distribution vector p.
  89. //
  90.   p=abs(h)/mysum(abs(h));
  91. //
  92. // Now we're ready to take some `walks'.
  93. //
  94.   total = 0;
  95.   for(i in 1:niters) {
  96. //
  97. // Generate an initial state for the particle, and then call randwalk.
  98. // Randwalk will keep track of the statistics we need for each walk.
  99. //
  100.     j = initstate(p);
  101. //  printf("Walk no. %d...",i);
  102.     eta = randwalk(A, P, f, j, h, p, chlen);
  103. //  printf("done\n");
  104.     total = total + eta;
  105.   }
  106.   return total/niters;
  107. }
  108.  
  109. //
  110. // End of file `estimate.r'
  111. //
  112.  
  113. mysum = function(X) 
  114. {
  115.   local(tmp, i, j);
  116.  
  117.   //
  118.   // Treat vectors, scalars and matrices differently.  If X is a scalar,
  119.   // do nothing...
  120.  
  121.   if (X.n == 1) { return X; }
  122.  
  123.   //
  124.   // if X is a matrix, return the vector of row-sums, except in the special
  125.   // case that X.nc==1 or X.rc==1; in that case, return the scalar sum.
  126.   //
  127.  
  128.   if (X.n  > 1)
  129.   {
  130.     if (X.nc == 1) {
  131.       tmp = 0;
  132.       for(i in 1:X.nr) {
  133.         tmp = tmp + X[i;1];
  134.       }
  135.       return tmp;
  136.     }
  137.     if (X.nr == 1) {
  138.       tmp = 0;
  139.       for(i in 1:X.nc) {
  140.         tmp = tmp + X[1;i];
  141.       }
  142.       return tmp;
  143.     }
  144.     tmp = zeros(X.nr);
  145.     for(i in 1:X.nr) {
  146.       for(j in 1:X.nc) {
  147.         tmp[i] = tmp[i] + X[i;j];
  148.       }
  149.     }
  150.   return tmp;
  151.   }
  152. }
  153.